class Gibbs_Thurst
{
	// Hyper parameter to be added
	public:
		Gibbs_Thurst(imat X);
		mat Z_update(mat m, double sigma, mat z);
		mat m_update(double sigma, mat Z);	
		double sigma_update(mat m, mat Z);	
		vector<mat> operator()(int M);
		mat InitZ(void);
	private:
		imat _X;
		int _n;
		int _p;
		mat _B;
		double _alpha;
		double _beta;
		
		rtnormNC *rt2;
		rtNC *rt;
		RandomG::Random<boost::random::normal_distribution<> > *rnorm;
};

 Gibbs_Thurst::Gibbs_Thurst(imat X)
{
	_X=X;
	_p=X.n_cols;
	_n=X.n_rows;
	rt=new rtNC();
	mat B(_p,_p);
	rt2=new rtnormNC();
	B.fill(1);
	_B=diagmat(B);
	_alpha=3;
	_beta=1;
	boost::random::normal_distribution<> N2;
	rnorm= new RandomG::Random<boost::random::normal_distribution<> >(N2);
}


mat Gibbs_Thurst::Z_update(mat m, double sigma,mat z)
{
	
	for(int j=0;j<_n;j++)
	{
		for(int i=0;i<_p;i++)
		{
			int ki1=-1;
			int ki_1=-1;
			int ki=_X(j,i);
			double t1=pINF;
			double t_1=-pINF;
			for(int l=0;l<_p;l++)
			{
				if(_X(j,l)==ki-1)
				{
					ki1=l;
				}
				if(_X(j,l)==ki+1)
				{
					ki_1=l;
				}
			}
					
			if(ki1!=-1) 
			{
				t1=z(j,ki1);
			}
			if(ki_1!=-1) 
			{
				t_1=z(j,ki_1);
			}
		//	cout << "t1 " << t_1 <<" " << t1 << " sif " << sigma(0,i)<< " m "<< m(0,i) <<  "\n";	
			if(t1==pINF)
			{
				z(j,i)=(*rt2)(m(0,i),t_1,pINF,sqrt(sigma));
			
			}else if(t_1==-pINF){
			
				z(j,i)=(*rt2)(m(0,i),-pINF,t1,sqrt(sigma));
			}else{
			
			 	//cout << "t1 " << t_1 <<" " << t1 << "\n";	
				z(j,i)=(*rt)(t_1,t1,m(0,i),sqrt(sigma));
				//z(j,i)=(*rt2)(m(0,i),t_1,t1,sigma(0,i));
			}
		//	cout << "z " << z(j,i) << " m " << m(0,i) << " s " << sigma(0,i);
		//	cout << " " << t_1 << " " << t1 <<" "<< ki <<" "<< z(j,i)<< "\n";
			
		}
	//	cout << "\n";
	}	
	return z;
}

mat Gibbs_Thurst::m_update(double sigma, mat Z)
{
	mat zb=sum(Z,0);	
	
	mat S(_p,1);
	mat me(_p,1);
	mat m(_p,1);
	for(int i=0;i<_p;i++)
	{
		S(i,0)=sigma/((_n+sigma/_B(i,i)));
		me(i,0)=zb(0,i)/(_n+sigma/_B(i,i));
		m(i,0)=(*rnorm)()*sqrt(S(i,0))+me(i,0);	
	}
	//cout << "m " << m.t();
	return m.t();
}
		
double Gibbs_Thurst::sigma_update(mat m, mat Z)
{
	double s=0;
	double a=_alpha+_n*_p/2;
	double b=_beta;
	for(int j=0;j<_p;j++)
	{
		for(int i=0;i<_n;i++)
		{
			b+=0.5*(Z(i,j)-m(0,j))*(Z(i,j)-m(0,j));
			//cout << "b " <<  b << "Z " << Z(i,j) << "m " << m(0,j); 
			
		}
		//cout << "\n";
	}
		cout << "a " << a << "b " << (double)1/b;
		boost::random::gamma_distribution<> Ni(a,(double)1/b);
		RandomG::Random<boost::random::gamma_distribution<> > rgamma(Ni);
		s=(double)1/rgamma();
	//cout << s;
	return s;
}

vector<mat> Gibbs_Thurst::operator()(int M){

	mat m(M,_p);
	m.fill(0);
	mat sig(M,1);
	sig.fill(1);
	mat z=InitZ();
	cout << z;
	for(int i=1;i<M;i++)
	{
		z=Z_update(m(i-1,span::all),sig(i-1,0),z);	
		m(i,span::all)=m_update(sig(i-1,0),z);
		sig(i,0)=sigma_update(m(i,span::all),z);
		cout << m(i,span::all) << ".";
	}	
	vector<mat> vec;
	vec.push_back(m);
	vec.push_back(sig);
	return vec;	
}

mat Gibbs_Thurst::InitZ(void)
{
	mat z(_n,_p);
	for(int i=0;i<_n;i++)
	{
		for(int j=0;j<_p;j++)
		{
			z(i,j)=(-(double)_X(i,j)+_p/2)/_p;
		}

	}
	return z;
}
